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Abstract. 


This reoort presents an Image method algorithm for the 
derivation of point sources of elastostatics in multi-layered medi 
assuming the infinite space point source is known. Specific cases 
have been worked out and shown to coincide with well known 
solutions in the literature. 
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Introduction 

The importance of point sources (Green’s functions are point 
sources) for some given (linear) governing differential equations 
and boundary conditions lies in two main reasons. First, any 
localized process when viewed from a sufficient distance can be 
modelled as some suitably chosen point sources. Second, Green's 
functions can be used to reframe the governing differential 
equations and boundary conditions in an integral equation form; the 
integral equation form can, for example, be used as the basis for 
numerically analyzing a large class of problems using the boundary 
element method. 

This paper presents an algorithm for the derivation of point 
sources of elastostatics in multi-layered med ia assuming the point 
source in infinite space is known. The method is similar to the 
image method that is familiar when deriving Green’s functions in 
plane 1 ayer ed med ia where there is only one unknown scalar field in 
the governing equations such as in temperature conduction, 
potential flow and electrostatics problems. The algorithm is then 
used to derive the Green’s functions for any point source in a 
region consisting of an elastic layer perfectly bonded to two 
elastic half spaces. 

Background 

There are many known Green’s functions for halfspace problems 
in elastostatics. Most of the known Green’s functions are 
specialized for a single halfspace having a stress free surface (a 
special case of bonded elastic halfspaces when one of the regions 
has zero rigidity). We will first briefly survey some of these 
known solutions with occasional comments on the method of 
derivation, then we will discuss some of the analytic methods used 
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to derive Green's functions In multilayered halfspaces, and mention 
a few known results for layered systems with some comments. 

The most used point source solutions are the point force, the 
dislocation and the nuclei of strain (or double couple) solutions. 
The point force solution for 2-D plane problems in a halfspace with 
a free surface (Mel Ian 1S32), 3-D problem, in a halfspace with a 

free surface (Mindlin 1936) and 3-D problem in bonded elastic 
halfspaces (Rongved 1955) are known. Rongved obtained the Green’s 
function through the use of the Papkovich-Neuber potentials and 
arguments from harmonic analysis; the resulting solution is in the 
form of the sum of a point force solution in infinite space and 
some point sources at the image point with respect to the interface 
plane. 

The screw dislocation (e.g. Rybicki 1971) and edge dislocation 
(e.g. Freund and Barnett ( 1976a, b) presented it in a convenient 
for,’Ti) for a .halfspace with a free surface and a general dislocation 
line intersecting a free surface (Yoffe 1961) are also known. The 
screw dislocation problem is obtained by the method of images 
(since there is only one field variable), while Freund and Barnett 
solved the edge dislocation problem through the use of complex 
analysis and the Mushkel ishvi 1 i potentials. 

There are six nucleii of strain sources. The solution to the 
first (double couple in a plane parallel to the free surface) was 
given by Steketee ( 1958), the re.maining five sources were given by 
Maruyama ( 1 964). Maruyama used image nudeii of strain sources to 
cancel the tangential component of the surface traction on the free 
surface. He then used the Boussinesq solution (in Galerkin vector 
representation) and the remaining normal tractions on the free 
surface in a Han kel /Fou r i er transformed space to obtain the rest of 
the fields afterwhich he t r ansf or.Tied the solution back to real 
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space. This procedure is highly specific to half space problems 

with a free surface and cannot be generalized to multiple layered 

♦» 

systems. We will next consider the known techniques to 
systematically treat sources in multiple layered media. 


A systematic formulation for the derivation of the fields due 
to 3-D sources in a layered halfspace was presented by Ben-Menahem 
and Singh (1968) and later refined by Singh (1970) and 
independently by Sato (1971). The formulation makes use of the 
analogue of Hansen's eigenvector expansion for electromagnetic 
problems (1935) applied to elastostatic and dynamic problems, 
combined with the Haskel 1 -Thompson transfer matrix technique and 
the Pekeris (1955) 'source condition’ at the level of the 
discontinuity. The formulation leads to a solution for the field 
variables of the form: 


O . . U . ae 

■i J , 1 


>d- 


Og'>< exp(p^kz) 


*»_ k^exp(6 kz ) 
n ^ 


J ( kr ) dk 
P 


( 1 ) 


where ; 

a , fi , 1 , t are constants dependent on the indices s 

s '^s q q 

and q 

J is Bessel's function of the first kind and 

P 

of order p 

r is the radial (cyl indrical ) coordinate 

z is the z cylindrical coordinate 

k is an integration variable 


In the above integrals, «.,•’£ snd £ may also depend on bot 
the elastic properties of the half space and the layer. 
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Sato and Matsu'ura (1973) and Jovanovich et al . (1974 a,b) 

made use of the above mentioned formulation to calculate surface 
deformations by numerically integrating the ensuing expressions. 
Such evaluations require special numerical techniques and 
significant amounts of computational effort. Furthermore, no field 
values were computed inside the layer or halfspace. To accomplish 
such evaluations requires significantly more effort (both in 
further algebraic manipulations and in computations and special 
numerical treatments); this observation 1s especially true for 
field points close to the source point. The importance of having 
the field variables being available everywhere occurs when a 
boundary element/integral equation formulation for processes 
occuring in a region consisting of such multi-layered media is 
r equ i red . 

A formulation in the same spirit as the Ben-Menahem and Singh 
formulation for multi-layered 2-D problems using the Airy stress 
function was presented by Singh and Gang (1985). Although the 
original formulation is applicable to 3-D, 2-D and antiplane 
problems, the specialized 2-D formulation is less complex. 

Simple'" but more specialized point source solutions in 
multi-layered media are available. For example, Rybicki (1971) 
presented the solution to a screw dislocation (can be specialized 
from his expressions) in a region consisting of an elastic layer 
with a free surface and perfectly bonded to an elastic halfspace. 
Rybicki used the method of images to derive his solution. Rundle 
and Jackson (1977) presented an approximate solution for a 3-D 
double couple source parallel to the free surface in a region 
consisting of an elastic layer perfectly bonded to an elastic 
halfspace. The approximate solution was obtained by using 
Steketee's 1958 solution combined with the use of an image method 
(using Rybicki 's technique) on the "antiplane" part of the point 
source being considered. 
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Rundle and Jackson compared thei'r solution for the surface 
deformations to that obtained by the direct integration of the 
improper integrals of Jovannovich et al. and found errors up to 14% 
while varying the thickness of the plate and the location of the 
source; however, the errors are less than 5% for the parameteric 
values they need in their specific application. We note that they 
take the rigidity contrast between the halfspace and elastic layer 
to be 10 to 1, a choice that would favor their a pproximat i on . 

Rundle and Jackson also obtained the (approximate) response of 
the point source they considered with a viscoelastic instead of an 
elastic halfspace through the use of the correspondence principle. 
The direct application of the correspondence principle to their 
elastic solution is possible because the material parameters are 
kept seperate from the geometric coordinates (in the form of a sum 
of material parameters multiplied by a function of the 
coordinates). This procedure is not directly applicable to the 
solution presented in the form of improper integrals (1) because 
the material parameters and the geometric coordinates are 
intermixed in the denominator of the integrand. 

Finally, we mention the existence of an image method for 
perfectly bonded elastic halfspaces in terms of the 

Papkovich-Neuber potentials (Aderogba 1977). Aderogba presents the 
algorithm for obtaining the four image potentials which involves 
multiple integrations with respect to the coordinate per pend i cu 1 a r 
to the interface plane and differentiation with respect to all 
three coordinates. The algorithm we present in this paper involves 
3 potentials only, and only differentiation of the potentials with 
respect to the coordinate perpendicular to the interface plane (as 
well as multiplication by scalars) is required. This distinction is 
especially important while contemplating the repeated use of the 
image algorithm to obtain the fields due to point sources in 
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regions consisting oH' an elastic layer perfectly bonded to two 
elastic halfspaces. 


Preliminary Cons iderat i ons 

The image method we consider here is dependent upon expressir. 
the displacements in terms of potentials. The specific potentials 
we use are the analogue to Hansen-’s potentials for el astostat ics 
and dynamics. Unlike Ben-Menahem and Singh (1968) we do not expand 
the potentials in terms of eigenfunctions, instead the algorithm 
operates directly on the potentials. Note however, that the 
derivation of the algorithm makes use of the eigenfunction 
expansion technicue (see Appendix 3). 


Specifically, we express the displacement field in terms of 
the Hansen potentials ¥* , and in the following way: 


= N(h,y.,) + F(a,h,l»2) + M(h,1»3) 



N(h,yp = VIP ^ (x , y , z-h) 

F(6,h,TP^) = 4 2-e^*J-1P,,(x,y,z-h) 

4 , Xi. O Z ^ 

- vY^(x,y .z-h) 

- 2-6*(z-h) v^¥’,,(x,y ,z-h) 

O /i, c. 


(2 


M(h,y»3) = V X 


le^Y ^{x,y 
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where : 


V is the gradient operator 

V X is the curl operator 


1 2 3 


A + 3fj 

is the Larne constant 


6 = 

A 

fj is the shear modulus 

h is a scalar for shifting the z-coordinate 


Note that the potentials and have to be harmonic in 

order for fi, £ and M to satisfy equilibrium. The Cartesian 
components for the displacements, strains and stresses are given in 
Append i x 1 . 


In order for these potentials to be useful for our purpose, we 
describe how to obtain these potentials given an elastic field 
satisfying equilibrium. We note the following: 


V • N = V • £1 = 0 


V X N = 0 


V • £ = 2 • ( 1 -a ) 


d^TP. 


dz‘ 


(3) 


- 

V X £ = 2* ( 1+e) i»e - 

dydz ^ 


2 . (1+e) 1-e 

dxdz y 


t X M 


"■’”3 - "•’^3 - "'’'3 

» e + — • e + , 


dxdz 


^ aydz y dz* ^ 


Therefore, if we have a given displacement field iJ, we 
calculate the following; 
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and 


7 - u = 2 • (1-5) 1 

az* 


a*TP 


o 


( V X u ) • e = 

" " az* 

From the above relations we find that: 


Y 


Y 


2 


3 




2 


V • u 
(1-6) 


+ 


z 



dz[(v 


X u ) • e 

“ z 


F2(x,y) + G2(x,y) 

+ z • Fg(x,y) + G 2 (x,y) 


(4) 


(5) 


( 6 ) 


The F.’s and G.'s are chosen such that Y- and Y^ are harmonic I 
11 2 3 I 

in the required region. Note, for the image method we should choose 
all the s ingu 1 ar i t i es of the potentials to occur in the region 
where the source occurs. This is made clearer in appendix 5 when we 
consider examples of the use of the algorithm. Finally, once Y^ and 
Y^ are determined, whatever remains is ascribed to Y^. If the given 
displacement field does satisfy equilibrium, the field should be 
expressible in terms of these three potentials (see Ben-Menahem and 
Singh 1968, and Morse and Feshbach 1953). 

The Hansen potentials for a point force, a line force 
perpendicular to the z-direction, and an edge dislocation 
per pendicu 1 ar to the z-direction are given in Appendix 2. We will 
show later (as is already known) that there is no need to obtain 
the potentials for a purely antiplane deformation field, since the 
image field involves a field of a si mi liar nature as the source. 
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The Algorithm 


Now we describe the algorithm and the notation associated with 
it; consider two elastic halfspaces perfectly bonded along an 
interface plane at z=0 (see figure 1). The material properties of 
region 1 are described by jj and 6 , and of region 2 by fj and 6^. 
Next we define the following: 

y(x,y,z) s T*(x,y,-z) 

^ 2 Hz/Hi ( 7 ) 

» = (6i + T) u = (Sx + n 

“ ■ (S^+y) ° • (ir-Cj + l ) 


Note that if Y is harmonic then IP is also harmon i c and hence 
can be used as a Hansen potential for £J., £ and M. 

The algorithm states that if we have the representation for a 
point source in infinite space of elastic constants similar to 
those of region 1 at the location x=y=0 and z=h described by the 
displacement field: 


0 _ 0,, . ^0 

y. = u (h,6j ,]Pp,E.. ) 


(S) 


then the displacement fields in regions 1 and 2 for a similiar 
point source in region 1 at x=y=0 and z=h are given by: 


1 

y. 



+ y. ( - h , 6 




u 


O 

c. 




( 9 ) 


where : 


Lr = BR(-h,a,b,6t).E.° 


Rl(0-Ll 




2 0 


gp^(-h,a ,b,6i ) = 


-26i(1-a)h*|- +(1-b) - 46?(1-a)h^-^ 

dz^ 


+(1-a) 


+26^ ( 1-a)h 


T^(h,a,b,6,,6i) B 



V . 3 1 

+ a 

-2(6jb-C»a)h-.^ 

0 

+ b 






We note that : 


= BR(-h,a,b,6,)-E.° = fp(-h,a,b,«»)-f.R 


a , b,6i ) = B|^( + h,a , b,6j ) 


Note that the ' are simple mu 1 1 i pi i cat i ves of The case 

when £.o 0 corresponds to the purely antn-plane problem, and thus 

R 

the algorithm reduces to the scalar image method for that case. 
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The derivation of the above algorithm is given in Appendix 3, 
an analytic check of the algorithm (making sure the displacement 
and traction interface conditions are satisfied) using the 
Cartesian components is shown in Appendix 4, and finally some 
sample known solutions are rederived in appendix 5; namely screw 
dislocation in a half space with a free surface, the Boussinesq and 
Cerruti point force normal and tangential (respectively) to a free 
surface, Flamant's line force normal to a free surface, a line 
force tangential to a free surface, and finally Mindlin's solution 
of a point force interior to a halfspace. 

In anticipation of applying the above algorithm to the 
derivation of point sources for a region consisting of an elastic 
layer perfectly bonded to two elastic halfspaces, we consider the 
effect of shifting the interface plane from z=0 to z=H on the form 
of the terms in the algorithm. 

Assume the interface is at z=H. Define a new coordinate 
z’ = z-H. If z = h is the location of the source point (h > H) then 
z’ = h-H is the location of the source point in terms of the new 
coordinate, and z' = H-h (z = 2H-h) is the location of the image of 
the source point with respect to the interface. In terms of z’, the 
algorithm is applicable as shown above with = h-H . We 

then reexpress z' in terms of z. Therefore, for the case when z=H 
is the interface plane we get: 

given: ii. = u (h,6j,2.f^,£,|^) 

then: + m( 2H-h , 6 , ^ ) 


= y.(h,C*,2.p,1P.L) 
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where ; 


2^-r = RpXH-h,a,b,6.)*£.° 








I^(-H + h.a,b.6j .«! )*]1 r 


2 0 
y>‘^ = T (y) -IP 
-L ^ -L 


( 12 ) 


Point sources nn a region consisting of an elastic layer perfectly 
bonded to two elastic halfspaces. 

In this section, we consider the method of derivation of the 
displacement field due to a point source in a region consisting of 
an elastic plate (0<z<H) perfectly bonded to two elastic halfspaces 
(see figure 2). The location of the point source is allowed to be 
either in the plate or in one of the halfspaces. The case of the 
source being in the elastic plate and the case of the source being 
in one of the halfspaces have to be treated separately. The elastic 
parameters used to characterize each region are chosen to be and 
6^., where 6^. is defined as 6^. = (h+fj) / (A+dfj) . 

In order to simplify the presentation , we define the following 
terms (this notation is suggested from private notes by Rice 1985 
of an outline of using the image method combined with the 
Papkovich-Neuber potentials to solve the same problem, although a 
detailed description of the implementation is not performed); 


+ 

H = 


^3 



+ 

6 



= 6 


. ± . 
= H 


6 


o 



16 


a* = (6 + 1 ) /(6 ^ +7* ) b* = 

(6^+1)/(^*6*+1) 

R*(h) s ,a*,b*) 


= 


i ✓ i 1 i i 

± 


± , 

T*(h) s I^(h,a=",b=",6 ,6^) 

II 


) 













a 5 + ) 

b 

“ 


* * 




ipCh) s Rp(h,«2,a ,b ) 

it 

“ 

Rl(- ) 

Xf X< X: 



* 

T^(h) s T^(h,a ,b 

II 

H 



Case I: source 1s in region 1 (the elastic plate'), and h<H 


(13) 


Referring to figure 3, a point source located in an elastic 
plate requires additional fields (derived using the algorithm) to 
conform with the boundary conditions on the upper interface. These 
additional generated fields have to conform with boundary 
conditions (derived using the algorithm at a shifted interface of 
location z=H ) on the 1 ower interface; hence each field generated by 
the algorithm requires a further image. The same argument applies 
when we start satisfying boundary conditions at the lower interface 
first. Hence, we can deduce the following: 


Given a point source in region 1 (h<H) in the form: 


= U(h,6 




) 


u = 


u" + u(-h,6^ (0)) 


m^1 


oo 

+ \ u(2mH + h,6 ,El"'^(ni),f, "^(m)) 

m^1 


+ u(2H-h,6^ ,f.R”(0),£.L (0)) 


+ } y.(2(m+1 )H-h,6^ ,]P^ (nn),r.^ (m)) 
m-1 


+ \ u(-2mH + h,6^ (m),]P|^'^ (m)) 

m^1 


2 

u = 


m^1 


u(- 2 (m- 1 )H-h, 62 ’E.R "’(rn)^£-L 


+ u(h,«2'2.R (O)-tL (0)) 


m’ 


00 

I 


y.(-2mH + h,6 1P^ (m),nf (m)) 


where : 


= SR(-h).£° 




2-r = Rp((2m-1)H + h)-f.^‘^(m-1) 


zl ^ ( m ) = 


" BR(-2mH-h)«E^ ■'(m) 


zl'^'Crn) = 


zl (0) = Bp(H-h).y° 


zl (0) 


^r"^ (m) = ER(-2fT>H + h)-L^ (m-1) 


IlI^ (m) 


B^S., 


flp (m) = BR((2m+1)H-h)*f.^'" (m) 


,1- 


n, (m) 


id (0) = Io(-H+h)-H° 


z; (0) 


ckcii oai ecu 
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S.R (m) = Tp^(-(2m+1)H + h).E.J'" (m) (m) = (m) 

£R”'^(m) = Ip(-(2m-1)H-h).Ep"'(m-1) ^ ^ ^ 


From the above recursive relations for the "image" potentials 
we would like to obtain their direct relation to the infinite space 


potentials ]f ^ . and ]f ^ . This can be done by induction and the final 
results are: 




m 

= I I fBp(-(2k+1)H+h) •Bp(-2kH+h)l -RlC-H+h)-?® 
[ k=1 I- ^ ^ Jkj ^ 

^(0) = e[-E° 
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r ^ 


- 


1 i 



r: 


R, ♦ R. 



1 1 

k=1 

= L =L 

k 



( 17 ) 


Once the potentials are known in terms of the 
2 

potentials, the IP potentials are obtained by one further matrix 
operation as given previously (16), 

Note that whereas the and matrices consist of one scalar 
per matrix, the R_ and matices are 4x4 matrix operators 

n 0 n H 0 2 

involving — and ^2 operators as well as constants. 

One subtle point when deriving any given point source for this 
case (i.e. when the point source is in the elastic layer) is that 
the potentials have to have all their singularities in the region 
z>0 for the series of "image" potentials generated with the first 
reflection being with respect to the upper interface, and the 
potentials have to have all their singularities in the region z<H 
for the series of "image" potentials generated with the first 
reflection being with respect to the lower interface. These 
conditions are imposed on the choice of the potentials defining the 
original source, in order not to introduce any further 
singularities in any given region thru the use of the algorithm. 
This condition can be satisfied due to the flexibility in choosing 
the potentials. 

Case II: source is in region 2 Tthe lower halfsoacel and h>H 

Referring to figure 4, a point source located in an elastic 
half-space requires additional fields (derived using the algorithm) 
to conform with the boundary conditions on the first upper 
interface (i.e. the elastic plate and lower halfspace interface). 
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i 

Fields in the elastic elate are thus generated, and these have to 
conform with boundary conditions on the upper interface (z=0) via 
the algorithm. From thereon each new field generated influencing 
the elastic plate have to have a further "image". However, there is 
only one "series" of sources generated by this process since the 
I first elastic field generated in the elastic plate already 

satisfies the boundary conditions with respect to the lower 
; halfspace. Hence, we can deduce the following: 

Given a point source in region 2 (h>H) in the form: 
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2- * 0 ^ 2- * 0 
zl (0) = BR(H-h).tR (0) = Rl-^l 

3PR"’^(m) = I~(-(2m-1)H-h)-V^'^(m-1) " X[ '“l"^ ^ ^ 

( 20 ) 

From the above recursive relations for the ’’image" potentials 
we would like to obtain their direct relation to the infinite space 

potentials IP ^ . and . This can be done by induction and the final 
results are: 


r.^^'CO) = gR(h)*I^(-h+H).lP° 




m 


k=1 


R^(2kH+h) •R^C(2k-1 )H+h) 
— ' K — ’ K 


— 1 4- -h * ^0 


R^(h)-I*(-h+H).]P° 






E;-"(m) = 


m 

1 

mr 

[r^.r- 

II 

11 

11 



k=- 

m 

n 

1 1 

k=- 

[[ 


+ * —0 

Sl-Il-El 


Bj^(-(2k-1)H-h)-Rj^(-(2k-2)H-h) 


TR(h-H).r-^ 


L-ib 


Jk 


= L -L 


(21 ) 


Again^ once the IP ^ potentials are known 1n terms of the 
2 

potentials, the TP potentials are obtained by one further matrix 
operation as given previously (20). 

When deriving any given point source for this case (i.e. when 
the point source is in the lower halfspace) the potentials have to 
have all their singul ar i t ies in the region z>H, in order for the 
reflected images not to introduce any further singularities within 
any given region. 


22 


Conclusions and further recommendations 

A vector image method has been presented, for elastic problems 
with planar layering. An algorithm has been presented on how to 
derive point source solutions for two bonded elastic halfspaces, 
and then extended to the case of an elastic layer perfectly bonded 
to two halfspaces. Specific cases have been worked out and shown to 
coincide with well known solutions in the literature. A study of 
the number of image potentials required to obtain a good 
approximation for the layered medium problem has not been 
investigated, but is of considerable practical interest. 

A feature of this image method is that solutions are obtained 
in either closed form or as infinite series of relatively simple 
expressions. The form of the expressions allow the use of term by 
term application of the cor res pondence principle for obtaining the 
viscoelastic response of the solution of the corresponding point 
source in a viscoelastic layered medium. 

The method of deriving the algorithm suggests that an 
analogous algorithm can be obtained for spherical interface 
problems, and 2-D (but not 3-D) cylindrical interface problems in 
el astostat ics . This suggestion is supported by the existence of a 
scalar image method and Hansen potential representat ions for both 
these geometries. 

Finally, it would also be of interest to Investigate 
equivalent algorithms for other governing equations. For example, 
el astodynamics and poroel ast ic ity could be potential candidates for 
such an investigation. El astodynamic problems, in particular, do 
have Hansen potential representations that have been well 
established and used and could be investigated first without the 
considerable preliminary formulations that are needed for 
poroelastic problems. 
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Appendix 2: Sample potentials for some point sources 


Point Force: 

The displacement field due to a point force at the origin can 
be written as : 


u . = 




) . + 6 • p, • X. 

1 r k 


1 

""i ? 

r 


where : 


6 S 


A+fJ _ 1 


1 


A+2fj ~ K 3-4u 


X + y 


+ z 


( 2 . 1 ) 


P 


i 


are the magnitudes of the point forces in the 
i ^ ^-d i recti on 


For the point force it can then be checked that the Hansen 
potentials are: 

"l =-§-•[ Pl-FII " ^ 2 -^ * P 3 --in(ptz) ] 

’’a = -T-[ -‘’l-FI? - P2-rlz * ] <^-2) 

>■3 = /!•[ ] 

where: p = 1 / [4n/i • ( 1 +6 ) ] 

Note that if the upper (lower) "sign" is chosen in one 
expression, the upper (lower) "signs" must be chosen throughout for 
all the potentials. Also note that taking r+z (r-z) in the 
expressions makes the potentials (but not necessarily the 


The displacement field due to a line force can be written (fo 
plain strain) as: 


1 I 

_ a 

^ -p^. • InC + 

* ■> 1 


for i,k = 1,3 

(2.3 

and 

= 0 

° P|< '‘k ’*1 2 


where : 

a : 

■ T+7JJ 

» * 

A +3// 





2^2 
= X + z 





Pi 

and Pj 

are the magnitude of 

the line forces 


For the line force it can be checked that the Hansen 
potentials are: 




1 87i/i6 


z*arctan(— ) - x*ln? + ( 1 +0 ) • x 


Po- 


Z’ln{ - z + x*arctan( 


|)] 




8 n /^6 


■p, • z*arctan(— ) - x*lnt + (1+6)*x 


~ p 


5 • [ z • Int - 


z + x*arctan(— ) 
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Dislocations parallel to the z =0 plane: 

The displacement due to a dislocation along the y-axis (plane 
Strain) can be written as: 


u . 

1 


^•arctan(£) - | 




and 


u, = 0 


for i , k , n = 1,3 


(2.5) 


where ; 


ik 


+1 for i = 1 , k = 3 

-1 for i = 3, k = 1 

0 otherwise 


d, and d.^ are the slip magnitude of the dislocations 

I kJ 


We note that the terms in the second brackets expressing the 
displacements are of the form of line force expressions with 
equivalent magnitudes of and thus their Hansen potentials 

are already known. The Hansen potentials for the terms in the first 
bracket can be shown to be: 


^Srac^et _ 1 

1 " 2^’ 



z*ln( - z + x*arctan(^) 
z*arctan(— ) - x»lnC + x 


.bracket 


= 0 


^Irac^et 

3 


= 0 


( 2 . 6 ) 
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Append -ix 3: Derivation of the image method algorithm 

The eigenfunction expansion method for elasticity problems in 
layered media was first formulated by Ben-Menahem and Singh 1968. 

We have used this method to derive the algorithm discussed in this 
paper. The notation (as far as possible) is the same as in the 1968 
reference paper, although some new temporary terms have been 
defined in order to simplify the algebra for this specific 
implementation . 

Any elastic displacement field satisfying the equilibrium 
equations : 


v*u + (1 +A/a/)*w*u = 0. 
can be written as the sum of fi, F. and M: 


N = 
F = 

fl = 







2-6 


z 



(3.1) 


(3.2) 


where; 6 = ( A +/^ ) / ( A + 3/i ) 

V*1P = 

Using the method of the separation of variables in cylindrical 
coordinates on the potentials T , TP ^ and T ^ in the form: 

I W 


We get: 
where : 


r> = R(r)-F(e)*Z(z) 

T = exp(±kz)*J ( kr ) • exp(±im0 ) 
m 

IS Bessel's function of the first kind order 


(3.3) 
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Reexpressing TP., TP and TP„ in the above form and carrying out 

the V and v x and operations, we get: 

o z 


m^1 


+ a -N + + B -F + + c -m 

mm mm mm mm mm mm 


dk 


(3.4) 


where : 


A*, B* and C* are constants dependent on ’m’ only, 
mm m 


^ = exp( ±kz ) • 

rri 


[ 


+ B 
m ~m 


F = exp(±kz)- (±1-2Skz)-P - (1±26kz)- 

~m m 


B 1 

-m J 


and : 


M = exp ( ±kz ) • Q_ 

“m ~m 


P = e *J ( kr ) • exp(im® ) 
““m z m 


R =: (e » ■ + e^’—^ — *4^)3 ( kr ) • exp( ime ) 

•“^m ^ r TRr e kr de^m'' ^ 


C = Ce •— - e«*T^)3 ( kr ) • exp( imfi ) 

^m ^ r kr d9 9 dkr' m' ' 

In the above expressions for P , B and C there is the 

m m m 

implicit understanding that we can consider either the real or 
imaginary components of the expressions seperately. 


(3.5) 


From the above expressions for the displacements, we can find 
the expressions for the tractions at a plane z=constant, and we 
rewrite the above as: 


where 
and : 

y = A^*exp(kz) + A *exp(-kz) 
m m m 

+ • ( - 1 -26kz ) • exp( kz ) + • (- 1 +26kz ) ♦ exp( -kz ) 

z = C^*exp(kz) + C 'exp(-kz) 
m m ' ' m 

Xfn = A^*/i*exp(kz) + A^*/i*exp(-kz) 

+ B^ •f/6 • ( 1 -2kz ) • exp( kz ) + B^ */i6 • ( 1 +2kz ) • exp(-kz ) 

'^m “ A^*^*exp(kz) - A^*^* exp(-kz ) 

+ B^ • /L/6 • (-1 -2kz ) • exp( kz ) + B^-/i6 • ( 1 -2kz ) • exp(-kz ) 

“ cl^^^-expCkz) - exp(-kz ) 
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R ' 

We notice that that u components are uncoupled from the u'" in 

m — m 

the sense that the A*, coefficients do not affect the u*" 

mm — m 

i R 

components and the C coefficients do not affect the u components. 

rn m 

R L 

We therefore treat the u and the u components seperately when 

m “m 

analyzing a specific problem in terms of the Hansen potentials. 

Now we consider the specific geometry shown in figure 3.1. The 
region consists of two elastic materials seperated by a planar 
interface. The material elasticity parameters used to characterize 
the regions are taken to be fj ^ , 6^ and f/j , 6a . A point source 
exists at the position z=-h. We are required to find the 
displacement fields for region 1 (z<0) and for region 2 (z>0) under 
the influence of the point source, such that the displacements and 
the tractions are continuous across the interface plane (z=0). 

R L 

In what follows, we are manipulating u and u in equation 

® ~m m 

3.6 for a fixed 'm', but the 'm' subscript will be dropped. First, 
we express the displacement and traction (on a z-plane) for 
(z+h) > 0 (which includes z=0) of a point source of arbitrary 
nature (using the eigenfunction expansion method and expressing the 
m'th component in matrix form) in the following way: 


y(E) 


-1 

-1 - 26jk*(z+h) 


^0 

y 

TCP) 


+ 1 

+ 2k^^^ 

-1 + 26ik*(z+h) 

+ 2k/i,6i-[1+2k*(z + h)] 


^0 

K§) 

0 

-2kiJ^ 

+ 2kA<i6i • [l-2k* (z + h)] 



exp(-k I z+h I ) 


+ 1 

-kfj^ 


Cq • exp( -k I z+h I ) 


and we define: 


u(P) 

T(B) 



0 z = 0 



u(C) 

T(C) 


Now the elastic fields in region 1 are expressible as 


u(P) 

T(B) 



+ 1 

+1 - 26,k*z 

+ 1 

-1 - 26jk-z 

+ 2kiJi 

+2k^i6i ' ( 1-2kz) 

+ 2k/ij 

+ 2kfJi6^ • (-1-2kz) 


• exp( kz ) 


T(C) 


+ 1 

+ k^. 


terms due to the 
point source 
as given above 


exp( kz ) 


terms due to the 
point source 
as given above 


And the elastic fields in region 2 are expressible as: 
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u(P) 


-1 

1 

N 

« ■ 
»o 

i 

I 1 

1 

1 i 

i 


^2 

y(§) 


+ 1 

~ 1 + 26 2 kz 


"2 

K?) 


+ 2k/L/j 

+ 2 kiJ 2 ^z • ( 1 +2kz ) 



T(B) 

. J 

2 

-2kiJ^ 

+ 2k/Jj6j • ( 1-2kz) 



exp( -kz ) 


u(C) 


1(C) 


J 2 



exp( -kz ) 


(3.13) 


Applying the condition that u and T are to be continuous 
^ ^ *Tn m 

(for each m) along the interface plane z=0, we get: 


+ 1 

+ 1 



+ 1 

+ 1 



+ 1 

-1 


+ 1 


®i 

+ 2kfJ^ 

+ 2kfJ ^6 y 

- 2 kfJy, 

- 2 k|^ 26 j 


_^L 

+ 2kiJ, 

-2kfJy6y 

+ 2k^2 

- 2 k/i 2^2 


. ^2 . 


+ 1 

-1 


1 


[Ssl 

+ klJy 

+ k/y. 


. S . 


. ^6 . 


(3.14) 

Now we solve for and A^ , by inverting the 

4x4 and 2x2 system of equations. We obtain: 


a; 


+ 4 / 2 - 7 - 6 i 

^A/ 2““7 —6 ^ 




1 

■^A/ 2 -^'i 6 i & 2^^ i 

-^/ 2 + 76 i 52 + 5 j 

762 + 1 

762+1 

_^L 

A 

+ 76^62 1 

1 

-76162-61 

762+1 

_ 2 kifx__ 

,^62+1 

.lEEn. 

^2 . 


+ 7 + 6 j 

+ 7 +6 1 

1 

7+6 1 

7+6 1 

J 


'- fL ' 

1 

+ 7 

_! Hk _ 



. s . 

7 + 1 

-1 

I 

. 


. ^6 . 


( 3 . 15 ) 

where; 7 5 

A s 2 - (7+6i ) • (762 + 1 ) 


Express 1 ng the S’s 1n terms of A^, and C^, and simplifying 
the expressions we get: 


A 

B 

A 


+ 

1 

+ 

2 

2 


0 

1 - b 

1 - a 

26i • (1-a)*kh 

a 

26 1 a • kh 

1 

0 1 

b 

i 

. 


exp(-kh) 




exp(-kh) 


( 3 . 16 ) 


where : 


a s ( 6 j+ 1 )/( 7 + 6 j ) 


( 3 . 17 ) 


b = (6 I + 1 )/( 7 • 6 2 + 1 ) 


Therefore we find that the displacement field (for a given m) 
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in region 1 and region 2 can be written as: 


U(P) 


+ 1 

+ 1-2fij kz 


0 

1-b 


’ ^0 ' 

u(B) 

1 

+ 1 

-1-26ikz 


1-a 

26i • ( 1-a) • kh 




exp [k(z-h) ] 


+ 


source terms J 


• « 
u(C) 


+ 1 

• 

1-7 

TT7 

• 

^0 

. 

1 

- 


- 


- 


exp [k(z-h) 3 


+ 


source terms J 

(3.18) 


u(P) 


-1 

-1-2Cikz 

u(B) 



-1+26i kz 

■ 

2 

• 

. 


-fa 

26^3 • kh 



0 

+ b 


B~ 




0 J 


exp[-k(z+h)3 






2 


- 1 


— 

+ 1 

• 

TfT 

• 

o 

o 

• 

2 

- 




L J 


exp [-k(z+h) ] 


(3.19) 


The y.(£) terms are in a form from which we can deduce the 
algorithm, however, the y.(£) and u(B.) terms have to be further 
manipulated. We now try to express the ii(£) and iJ(S) terms in the 
following manner; 


u(P) 

u(B) 

■I 1 


+ 1 

+1-26ik*(z-h) 


a 

-f 1 

-1-20ik*(z-h) 


6^ 

a 

1 

+ 1 

+ 1-28j k- (z-h) 



+ 1 

-1 -26 i k* (z-h) 




exp[k(z-h)3 


a 

a z 


• exp[k(z-h)3 
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^ _a! 

+ 1 

.[<] 

•exp[k(z-h)] 

+ 

source terms 

dz^ 

1 

+ 1 

r 1 


1 r 

1 

(3. 


u(P) 

u(B) 


d 

dz 


dz‘ 


exp[-k(z+h)] 


-1 

-1-26ik- (z+h) 


A 

a_ 


-1+26ik* (z+h) 


b" 

a 

“1 

-1-26 jk* (z+h) 


i^l 

+ 1 ' 

-1+26jk' (z+h) 


1 

1 SI 
00 

« 

r 

1 r - 1 




exp [~k(z+h) ] 


exp [-k( z+h) ] 


(3.21) 


Not 1 ng that ; 

.n 

exp[k(z-h)] = k *exp[k(z-h)] 

dz'^ 

and 

.n 

exp [-k ( z + h ) ] = ( -k • exp [-k ( z + h) ] 

dz'^ 

We obta in: 

< = (l-b)-B- b; = (l-a).A- 

= -2Si • (1-a)-h-A~ + 45? • ( 1-a) • h- b" 

B^ = 26i • ( 1-a) • h- B~ A^ = -46? • ( 1 -a ) • h^ • B~ 


(3.22) 


(3.23) 
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and : 


A — a • A « 
a 0 


8 = b*B_ 

a 0 


(3.24) 


Aj^ = 2 • (6 2 b-6 j a ) • h • Bq 


B, = A = 0 
b c 


Since the above relations are true for each component of a 
potential, then they must be true for the whole potential and we 
get : 


if : 


= N(-h,1-°) + F(6»,-h,y»°) + liC-h,^®) 


Then : 


= u° + N(h, (1-b)-r>°) + |^N(h,-26,.(1-a)-h-F°) 

+ |jN(h,+46T (1-a).h-F2) 

2 

+ -^N(h.-46i .(1-a)-h^.y°) 


dz 

+ £(6, ,h, (1-a)F°) 

- M(h,4^.F°) 


(3.25) 


liC-h.aV’^) + |^E(-h,2-(6ab-6,a)-h-r»2) 
+ §7£(6 2 ,-h,b1P°) 


(3.26) 
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We obtain the algorithm given in the main body of this paper 
(two minor differences are: i) The statement of the algorithm in 

the paper considers region 1 to be at z>0 and hence z=+h instead of 
z=-h to be the location of the source point and ii) A formalism in 
terms of matrix operators is implemented in the main text). 
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Appendix 4; Analytic check of the image method algorithm 


Referring to figure 1 and using the notation defined in the 
main text, the image algorithm states that if: 


u = 


Then : 


N(h,TP°) + F(6i,h,1p2) + M(h,F°) 


u° + u\-h,«j 


y. = 


2 ,. _ ^2 _ 2 . 
u (h,6j 


where : 


e ; 

4 

4 


Bf^(h,a,b,6j ) -y® 
Ip^(h,a , b,6j ,6^ ) 

T,0’)-£° 


(4.1) 


We have already expressed the Cartesian components of the 
displacement and stress fields of any given Hansen potentials 
(Appendix 1). In this section, we will check whether the Cartesian 
components of the displacement and stress fields (in terms of 
Hansen’s potentials) that are generated by the image algorithm 
satisfy the conditions of displacement and traction continuity 
along the interface plane. 


In order to simplify the checking of the algorithm, we will 
consider the following 3 cases seperately: 


T) TP, = TP^: 


ii) TP^ = TP ; 


= 0 


m ) y, = 1* ; 


y, = yj = 0 


We will also note the following: 


a — r a ■ 


:> |_y> 

az 


z=0 


-Ir^ 

dz 


z=0 


(4.2 


l-y 

ax 


z = 0 


ax 


snd 


z = 0 


i-F = i-F 

z=0 2=0 


Also def i ne : 


= 


^F° 


^(x,y,z) a(x,y,z) 'z=0 


(4 . 3 


^F 2 


^0 


d(x,y,z) a(x,y,z) *z=0 


Case I : y = y " : y = r . = 0 

\ 4^ sj 


'z = 0 


®x‘ 26i.(1-a)-h*^^y + (l-a).(-i^y - 26-h-^^r) 


dxdz 


dx 


dxaz 


2 2 

+ e • fl— y + 26j • ( 1-a) • h'^^^y + (1-a)-(-|— F - 26»h- ,^, F ) 
y [ay * ' ^ ayaz ' ^ dy dydz 

2 . .2 


+ e 


l^y + 26 , . (l-a).h.-ilF + (l-a)-(i^F - 26 .h.-l^y) 


z dz 


dz 


dz 


(4.4 


= e 


z=^0 


[ *•1^’' ] ®y-[ ^•57’" ] * ^2-[ ^•57’’ ] 


(4 . 5 


We notice that the displacement field is continuous across z= 


(the interface plane). Next we consider the ± traction continuity 
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(equi 1 i bn’um) . 


-Zl 


= e 


2 = 0 


X 


2 3 

26, « (1-a) •2Afj «h t—W 


dxdz 


- 

+ (1-a)* (-2 a/, 6, .-_r - 4ju,6,.h- 


dxdz 


^>] 


+ e 


2 3 

f ■" 26 ,-( 1 -a). 2 iW,-h t—r 

*• ^ dydz 


2 3 

+ (1-a).(-2Af,6, - 4A/,6,-h i^F) 


dydz 

- r d^ 

+ e^*| + 26 ,.( 1 -a)* 2 /i,-h*-^y 


dydz 


dz 


dz 


2 3 

+ (1-a) • (2/i/,6, - 4)U,6, .h--i^r) J 


dz 


dz 


(4.6) 


-z2 


'z=0 




* %■( 


+ e 


[ ] 


(4.7) 


Noting that 


2 /^ 2 * [1 + (1-a)-6,] = 2^a-a 


s 1 nee 


(4.8) 




7 + —6] 


7 -f 6 


= 2^12 


6 , + 1 
If +6 . 


We find that the tractions across the interface plane (z=0) 


are conti’nuous. 


Case II: ^ y j £3-^-0 


= e 


' z==0 


z-0 


“x 


ax 


-y + 2C. h- 
< * dxdz 


(l-b)-i-y - 4 S^( 1 -a)•h^ — L— 

dxdz‘' 


_ 


ax 

26 ^ -h 


dxdz 




+ e 


2 

+ 26.h-j^r H- - 4»?.C1-a).h 


= -C? 


4 ^ _ 

26.-(l-a)-h-(-j^y - 28. -h 


aydz 




dydz* 


+ e 


2 3 

^ + (l-b)’l^y - 4«1 • ( 1-a) • h^-^y 

dz 


[ 7 f * TT 

*• dz 

2 3 

- 26 ^ • (i-a)-h- ( -i^y - 2ex-h--i^y) 

dz 


dz' 


(4 . 9 


= e 


[- 


a a 

2 • h* (626-618) + b*(“^ 1 P + 262 
^ * dxaz dx * 


axa 




[- 


+ e • -2 • h • (6 2 b-6 1 a ) 




3 yaz 

.2 


b*(-|-y + 2Ca-h--i— y) 
dy * dydz 


+ e • -2 • h • (6j b-6 1 a ) • 75-y + 


dz 


b-(j^y + 26* 




dz 


(4.10) 


We notice that the displacement field is continuous across z-0 

(the interface plane). Next we consider the traction continuity 
(equ i librium) . 






•y + 4 //j 



+ 2 A<i 


O-b) 



a ^ 
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+ e 






dxdz 




dxdz ' 


y - 4 Aii 6 i -h 


a^ ^ 


dxdz 


•TP) 


2 

+ • I + 4A^i-6^.h 


dydz 

46 ?-( 1 -a)* 2 f^,-h^' 


dydz 


y-y + 2A<i-(1-b) 


a 2 ^ 


aydz 


^y 


aydz ' 


26 ,.( 1 -a)*h.(- 2 A/i 6 , 1 -^ 

ayaz 


y - 4/ii«i • h' 




aydz' 

2 3 2 

2/^i6j*-^y + 4/ii .6i.h*-i^y -t- 2^1 •(i-b)--^y 

dz dz dz'^ 

2 2 ^ ^ 

- 467-( 1-a)-2Aix-h --^y 


dz 


- 26 x-( 1 -a)*h-( 2 A/i 6 x 


a 3 


4 /ixfi: 


az ' 




( 4 . 11 ) 


.z2 


= e 


'z = 0 


-2h • (62 b-6 xa ) • 2/^2 


axaz ' 


+ b- (- 2/^262 -;3^y + 4A/j6a*h 1-^y) 


+ e 


l- 


2 h- ( 62b-6 xa) • 


dxdz 


dydz' 


d^ .3 ■, 

+ + 4 ^ 26 *-h — 2_^y) 

^ dydz J 


•<- e 


-2h- (6 2b-6xa) • 2 a<2 

dz"^ 

2 3 

+ + 4 /J 2 C 2 *h--^y) 

dz"^ dz 


( 4 . 12 ) 


Noting that: 


+ (1-b)] = 2/ia-b-6* 


s 1 nee 




1626 I + i& 2 


1 + ir6. 


Sifi , ■*'®2 


1 + -»6j 


and 


4/i,[6, + 6t'(1-a)] = 4AfjC,-a 


(4.13) 


We find that the tractions across the interface plane (z-0) 
are continuous. 


Case III: y > 3 = y "; = 0 


= e 


z = 0 


+ e 


f l-y> . ;n.^? 1 

[ ay 1+7 ay J 
[ ax 1+7 ax J 


+ e 




( 4 . 14 ) 


= e 


'z==0 


X 


f — ■ 

[ 1+-y 


ay 


+ e * 

y TT 7 ax 


2 a 


+ e 


[° 


(4.15) 


We notice that the displacement field is continuous across z=0 

(the interface plane). Next we consider the traction continuity 
(equ i 1 ibr ium) . 
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+ e 


^ 3- 

<3xdz T+T dxdz 


+ e 


[°] 


(4.16) 


-z2 


= e 


z = 0 


+ e 


•[ 

[- 


/i* 


l+lf dydz 
2 


2 d 


1+7 dxdz 


[°] 


(4.17) 


Noting that 


1 - 


1-7 

1+7 


1^2 


1+7 


(4.18) 


We find that the tractions across the interface plane (z=0) 
are continuous. 

Therefore, all three cases check correctly, and the algorithm 


is correct . 
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Appendix 5: Derivation of some sample Green’s functions 

through the use of the image algorithm 

In this section, we consider displacement fields in Cartesian 
components for some sample point source problems on or in a 
halfspace with a free surface. The solutions that will be rederived 
are readily available (and established) in the literature and hence 
serve as an empirical check (see appendix 4 for an analytic check) 
of the algorithm. In addition, these specific examples help clarify 
details of the application of the algorithm. 

By considering a halfspace problem with a free surface, we 
obtain the following simplifications; 


-1 = 0 


call: 


then : 


Rp^(h,a,b,6i ) 


1-a = -1/5i 


6 s 




P .2 1 


-2h-|— 

dz 

-6 4- 46h • 7^ 


a 2 


-1/6 

+2h-|- 

dz 


. ' 

, 


1-b = -6j 




1 


(5.1) 

The following example problems will be considered: 

I. Screw dislocation (Antiplane problem). 

II. Line force (x=0, z=0) acting on a free surface. 

i) The line force is in the x direction. 

ii) The line force is in the z direction (Flamant's 
solution) . 

III. Point force acting on a free surface. 

i) The point force is in the x direction (Cerruti’s solution). 
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1i) The point force is in the z direction (Boussinesq 
sol ut i on ) . 

IV. Point force acting interior to the halfspace with a free 
surface (Mindlin's solution). 

i) The point force is in the x direction. 

ii) The point force is in the z direction. 

I . Screw dislocatio n (antiplane problem) 

For the antiplane problem, all that is required is to obtain 
the image potential with respect to the interface plane since the 
matrix is the identity operator. Also, we notice that getting 
the image of a given M type (see main text) displacement field is 
equal to the M displacement field of the image potential describing 

that field (i.e. M( h ) = fil(-h,E))» and hence we can directly 
operate on a given displacement field when using the algorithm for 
a purely antiplane problem. This corresponds to the scalar field 
image method for the antiplane case. 

As an example we consider the field due to a screw dislocation 
in the plane perpendicular to the x-z plane at location z=h and 
x=0. The field due to the dislocation in infinite space is: 

u^ = arctan [ (z-h)/x] (5.2) 

The image field will be which implies that the combined 
fields give: 


u^ = arctan [ (z-h)/x] - arctan [ (z+h)/x] (5.3) 

Of course this is but a simple application of the scalar image 
method . 
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1 1 . Line force (x=0.*z = 0) acting on a free surface. 

For this problem the potentials are given in appendix 2. We 
note that = 0 (i.e. there is no antiplane mode in a plane strain 
problem as is trivially known). 

We note that for the case when the point source is at the 
interface (i.e. h=0), we get a significant simplification in the 

matrix operator in the following way: 


Bp(0 ,a , b,6i ) 



0 

-6 


1/C 

» 

0 


(5.4) 


This means that for the displacement fields all we need to 
calculate are the following (e.g. see equations (2)): 

a a 1 [-0] fd d a^ 

dx’ dzj*[ 2j [ax' dz' axdz’ ^^2 

Now we perform the above dif f erentiations for the image 
potentials due to a line force. These potentials are linear 
combinations of the following functions: 



= ( z • arctan ( z/x ) - X'lnC + (1+6)*x)/2 = 
= (z*lnt - z + X • arctan (z/x) )/2 = 


2 2 2 

where: { = x + z 

We get: 



A 



C 


(5.5) 


-lnt/2 + 6/2 


[arctan(z/x)]/2 
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2 

= (z/e^)/2 

( 5 . 6 ) 

1) Case when the force 1s acting In the x direction: 

We get : 

0 

u. - u, 

+ [a/ ( 47T/75 ) ] • ^ 

e •[-6*1n?/2 - lnt/(26) - z^/t^ + constants ] 

X 

® “ ai^ctan ( z/x)/( 2C ) + xz/t^] j 

( 5 . 7 ) 

Noting the following identities; 

z^/e^ = 1 - x^/e^ 6/2-1/C26) = (1+6)/(A+^) 

(5.8) 

6/2+1/(2«)+1 = (1+6)^/(26) a = 2-6/(6+1) 

We find that the above solution coincides with that given in 
Love 1927 (article 151), except for a rigid body motion. 

i) Case when the force is acting in the z direction (Flamant’s 
solution) : 


= [arctan(z/x)]/2 
= -(z/C^)/2 

dz 


We get : 
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u = u 


- [a/(4Tt/j6 




e / [6 • arctan(z/x)/2 - arctan (z/x )/( 26 ) - xz/t^" ] 


+ e^-[6-Tnt/2 + lnt/(26) 


O O 

z'/r] 


(5.9) 


Noting the same identities mentioned above (5.8), we find that 
the above solution coincides with that given in Love 1927 (article 
151), except for a rigid body motion. 

III. Point force acting on a free surface. 

For this problem the potentials for the point source in 
infinite space are given in appendix 2. However, we have a choice 
of where to locate the singularities of the potentials. Since we do 
not want the image potentials to introduce any new sources inside 
the halfspace, we choose the infinite space potentials to have all 
their singularities in that ha If space. 

Again, if we are only interested in the case when the point 
force acts on the free surface (h=0), we get equation (5.4). This 
means that all we need to calculate are the following: 


a 

a 

d 



ax’ 

dy’ 

a7 

■ 

^2 


a a a a^ a^ a^l [-o] 
dx' ay' dz' axaz' ayaz' 


Now we perform the above differentiations for the image 
potentials due to a point force. These potentials are linear 
combinations of the following functions: 
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= x/[2-(r + z)3 
= -[1n(r+z)]/2 


(5.10) 


In addition, we have an antiplane potential for this case 
which Is a linear multiple of the following function: 


IPg = y/[2-(r + z)] 

and we will also have to calculate: 


(5.11) 


' a a 
ax' dy 


We get : 


= 1/(r+z) - xV[r.(r+z)^] j/2 


—IP = 
ay A 


■(xy)/[r- (r + z)*"] 


/2 


az’^A 




x/[r*(r+z)] /2 




dxaz A 


^ -1/[r*(r+z)] + x^/[r^* (r+z)] + x^/ [r ^ • ( r+z ) ^] j/2 


d ^ r 

-F. = (xy)/[r^- (r+z)] + ( xy )/ [ r ^ • ( r+z ) ^ ] 


ayaz A 


/2 




x/r ' 


/2 


(5.12) 


Ljf = r 

ax c [ 


•X/ [r • ( r+z ) ] /2 


axdz C 


■1^,. = x/r 


/2 


(5.13) 




~ ZZ \ 

dydz'^C [ 




4*^0 " f 1/^ 


Ijfg = I -(xy)/[P' (r+z)^] j/2 

^?g = [ l/(r+z) - y^/[r'(r+z)^] J/2 

■i ) Case when the force 1s acting 1n the x direction: 
We get : 


(5.14) 


0 

u = u 

+ [ 1/4ti)l/( 1 +5)] • 

e • 6*/[2*(r+z)] - (6 • x^ )/[2r • (r+z )^] 

2 2 

+ 1/[26-(r+z)] - X /[26r • (r+z) ] 

- z/[r*(r+z)] + (zx^)/[r^* (r+z)] 

+ (zx^)/[r^-(r+z)^)] 

+ (1+6)/(r+z) - ( 1+6) -y^/Cr • (r+z)^] 
+ e^- -(6-xy)/[2r • (r+z)^] 

- (xy)/[26r- (r+z)^] 

+ (xyz)/[r"^ • (r+z)] 

+ (xyz)/[r^- (r+z)^)] 

+ ( 1 +6 ) • (xy )/ [ r • ( r+z ) ^] 1 
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original page is 

^ POOR QUALfTY 


+ e 


-(6*x)/[2r- (r+z)] 
+ x/[26r- ( r+z ) ] 


+ ( xz 


<z )/r'^ j 


(5.15) 


Noting the following identities: 

(x^z)/[r'^* (r+z)] + (x^z )/ [r ^ • ( r+z ) ^] = x^/r^ - x^/ [r • ( r+z ) *■ ] 
z/[r*(r+z)] = 1/r - 1/(r+z) 

y'"/ [r • ( r+z ) “^ ] = 2/(r+z) -1/r - x^/ [ r • ( r+z ) ‘" ] 

(xyz )/[r‘^ • (r+z)] + ( xyz ) / [ r ^ ♦ ( r+z ) ^ ) ] = (xy)/r‘^ - (xy )/ [r • ( r + z ) ^] 
- 6/2 + 1 /( 26 ) = jLi’ (U6)/(A+JJ) 

(5.16) 

We find that the above result coincides with the published 
results (e.g. Love 1927, page 243). 

ii) Case when the force is acting in the z direction; 

We get: 


u = ii 


- [1/4n/u(1+6)] 


-(6«x)/[2r*(r+z)] + x/ [ 25 r«(r+z)] -(xz)/r‘ 


(C-y)/[2r- (r+z)] + y/ [26r • ( r+z ) ] - (yz)/r' 


(5.17) 


cerm’ng material 
t coincides with the 
published results (e.g. Love 1927, page 191). 

I V . Point force acting interior to a halfsoace with a free surface. 

The infinite space source potentials for this case are the 
same as for case III. However, since h^O in general, the matrix 
operator involved in the algorithm (5.1) requires us to further 
calculate in addition to the terms shown in case III, the following 
functions: 


£_1 \ w ^] 

dx' dy’ dzj ■ [ 1 J 


dz dx ' dy' 


fi_ £_ L-] fF°l 

.2 [dx' ay' azJ-[ 2j 


a a 

dxdz * dydz ' 



There are some repetition in the suggested functions to be 
calculated since the partial differentiation operations are 
commutative when the function is sufficiently smooth. Now we 
perform the additional required differentiations: 




3-x^/r^ 1/2 


. 2- A 
dz dy 


■3 • (xy)/r /2 


(5.18) 
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O 

-Cr. 


dr' 


-3-(xz)/r" /i 


a" - 


dr dx 


•3 • (xz )/r 


/< 


2. 


dr dy 


•3 • (yz )/r' 


/2 


[ 

dz '■ 


^ — [l/p -~3*2:/r' 1 / 4 : 


/ o 


t) Case when the force is acting in the x direction: 


(5.19) 


Def in ing ; 

000 O 

- X +y +(z + h)‘' 

We get after simplification and the use of identities similar 
to those given in (5.16), but with z replaced by (z+h) and ”r" 
replaced by "r' 2 " wherever they occur: 


U = 


0 


u 


+ [1/4h/l/( 1+6)] • 


e 


X 


+ e • 

y I 


L 

[6/2 + 1/(26) + 1 - 2-(1+6) + ( 1+6)]/(r2+z+h) 

+ [-6/2 - 1/(26) - 1 + (1+6)] .x^/[r2*(r2+z+h)^] 
+ [-^ + (1+6)]/r2 

+ x^/r^ + 26 hz*[ 1 /r 2 - S-x^/r^] 

[-6/2 -1/(26) - 1 - (1+6)] .(xy)/[r--(r^+z+h)^] 

+ (xy)/r“ + 26hz-[1/rJ - 3-x^/r^] J 


+ [-6/2 +1/(26)] -x/Cr^- (r^+z+h)] 

+ {z-h)'x/r^ - 66hz • X* (z+h)/r2 j 

(5.20) 

Noting t hat : 

6/ [411 A#- ( 1+6)] = 1/[16wa/- (1-1^)] 

1/6 = 3-4v (5.21) 

-6/2+1/(26) = 46 • ( 1-1^) • ( 1-2y ) 

We find that the above result (5.20) coincides with the 
solution first obtained by Mindlin (1936) and shown in Mura (1982). 

ii) Case when the force is acting in the z direction; 

After simplifications we get: 


u = u 


+ [ 1 /4ha/( 1 +6 ) ] 


+ e 


+ e 


+ [6/2 - 1/(26)] -x/Lr^* (r^+z+h)] 


3 5 1 

+ (z-h)’x/r2 + 66hz • X • ( z+h ) /r 2 1 


[6/2 -1/(26)] •y/[r2* (n^+z+h)] 


+ (z-hO-y/r^ + 66hz • y • (z + hO/r^l j 


[6/2 +1/(26)]/r. 


[(z+h)^ - 26hz]/r2 + 66hz • (z+h)^/r2 j 


(5.22) 


Noting the relations given in (5.21) and: 



60 


6/2+1/(26) = 6.[8*(1-v)^ - (3-4P)] (5.23) 

I We find that the above result (5.22) coincides with the 

I 

solution first obtained by Mindlin (1936) and shown in Mura (1982). 


Figure Captions 
figure 1 ; 

2 bonded elastic halfspaces with a point source at z=h. 
figure 2 : 

An elastic plate of thickness H perfectly bonded to two 
elastic halfspaces, with a point source in region 1 (the elastic 
plate). 

figure 3 : 

Location of the point source and image point sources for the 
elastic field in region 1 (left series of points) and region 2 

(right series of points) when the point source is located in region 

1 (the elastic plate). 

figure 4 : 

Location of the point source and image point sources for the 
elastic field in region 1 (left series of points) and region 2 

(right series of points) when the point source is located in region 

2 (the lower halfspace). 


figure 3.1: 

2 bonded elastic halfspaces with a point source at z--h. 



Figure 1 
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Figure 3.1 


